Phase behaviour of colloidal assemblies on 2D 
corrugated substrates 

Samir El Shawish 1 , Emmanuel Trizac 2 and Jure Dobnikar 3,4 

1 Reactor Engineering Division, Jozef Stefan Institute, Jamova 39, 1000 Ljubljana, 
Slovenia 

2 Universite Paris-Sud, Laboratoire de Physique Theorique et Modeles Statistiques 
(CNRS UMR 8626), 91405 Orsay Cedex (France) 

3 Department of Theoretical Physics, Jozef Stefan Institute, Jamova 39, 1000 
Ljubljana, Slovenia 

4 Department of Chemistry, University of Cambridge, Lensfield Road, CB2 1EW, 
Cambridge, UK 

E-mail: jd489@cam. ac . uk,trizac@lptms . u-psud. f r 

Abstract. We investigate - with Monte Carlo computer simulations - the phase 
behaviour of dimeric colloidal molecules on periodic substrates with square symmetry. 
The molecules are formed in a two-dimensional suspension of like charged colloids 
subject to periodic external confinement, which can be experimentally realized by 
optical methods. We study the evolution of positional and orientational order by 
varying the temperature across the melting transition. We propose and evaluate 
appropriate order parameters as well as the specific heat capacity and show that the 
decay of positional correlations belongs to a class of crossover transitions while the 
orientational melting is a second order phase transition. 



PACS numbers: 82.70.Dd, 82.20.Wt. 83.80.Hj, 79.60.Jv, 81.07.Bc, 82.70.Kj 

Submitted to: J. Phys.: Condens. Matter 
1. Introduction 



Understanding, creating and controlling patterned structures on nano- and microscale 
is an important objective of modern materials science pQ. Studying colloidal ordering at 
surfaces provides valuable insights into the underlying mechanisms. Complex surfaces 
with submicron periodicities can direct self-assembly of highly oriented responsive 
colloidal crystals [2] and biologically active substrates [21 0]. Charged objects are 
ubiquitous in colloidal domain and electrostatic interactions have been widely employed 
in colloidal assembly. Predominantly spherical colloidal particles have been studied, 
while comparatively little work has been devoted to the behaviour of anisotropic charged 
composite objects in a solution. The spherical shape is, however, more an exception than 
a rule in the colloidal realm, and from the fundamental point of view, understanding 
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the mechanisms directing the assembly of anisotropic charged colloids is at least equally 
important. Moreover, it enables design of novel ordered micro- and nanostructures. 
Experimentally, patterned substrates are typically created by either optical manipula- 
tion [5 J or surface deposition [6j [7J. Multiple charged spherical colloids can then be 
confined within the traps and if the number of colloids iV is an integer of the num- 
ber of traps N tT , at strong enough confinement, monodisperse clusters or "colloidal 
molecules" (in general anisotropic, at least for small n ) constituted of n = N/N tr col- 
loids are confined to each of the traps. The coupling between the anisotropy of the 
charged "molecules" and screening by microions in the solution is surprisingly nontriv- 
ial p2J ttH US] ■ The multipole moments generated by the molecular anisotropy are all 
screened in a similar fashion, so that inter-molecular interactions feature non trivial 
angular dependence affecting the symmetry of the ordering in a subtle way. Therefore, 
at large enough electrostatic interactions, orientationally ordered "colloidal molecular 
crystals" are formed. Crystals with remarkably rich variety of orientational ordering 
have been observed P El HQl HU [121 H31 HH US HH1 [IT]. In the case of like-charged 
repulsive colloids the molecular size is determined by the interplay between the elec- 
trostatic repulsion and the confinement forces, while the intermolecular electrostatic 
interactions determine the type of rotational ordering of the molecules. The ground 
state behaviour of such colloidal molecular crystals has been explored in several studies 
P2J [El [111 [15j [TFJ , however, in order to assess their thermal stability the complete phase 
behaviour needs to be explored. In [10] molecular dynamics simulations have been per- 
formed showing the generic phase behaviour of colloidal molecular crystals (n = 2, 3 and 
4) on square and triangular periodic substrates. A two-stage melting has been reported 
where first the orientational and then the positional order vanishes. Similarly, in an 
experiment [U [9] laser induced freezing and melting has been observed: upon increasing 
the confinement strength colloidal trimers confined to triangular lattice went through 
a transition from liquid to crystal and back to modulated liquid. Here, we focused on 
colloidal dimers (n = 2) on two-dimensional square periodic substrates and performed 
Parallel Tempering Monte Carlo simulations to study the melting process in detail. We 
defined appropriate order parameters describing the two-stage melting and character- 
ized the order-disorder transition for both, the orientational and positional correlations. 
It turns out that the orientational melting is a second order phase transition with di- 
verging specific heat at the transition temperature, while the decay of the positional 
order is not a phase transition but rather a gradual crossover-type transition. 

2. Model and methods 

We consider a system of like charged spherical colloids confined to a plane and immersed 
in a suspension of monovalent counterions and salt ions (the solvent is hereafter 
considered as a structure- less medium of constant dielectric permittivity). The microions 
neutralize the total charge and mediate repulsive screened Coulomb interactions among 
the colloids. We assume the pairwise additive Yukawa potentials and do not consider 
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Figure 1. Schematic representation of colloidal dimers in potential traps. The 
distance between the centers of the potential traps I is the unit of lenght. In strong 
enough confinement, exactly two colloids occupy each trap forming a colloidal dimer 
with a characteristic size 2d. The colloidal "molecules" are well defined, when 2c? <C I- 

many-body contributions to the electrostatic interactions because the geometry of the 
problem and the typical screening parameters explored are in the regime where pairwise 
additivity is an excellent assumption [18]. The colloids are subject to additional periodic 
external confinement, introduced by laser tweezers [H [10] creating a set of potential 
traps for the colloids. The number of traps iV tr = N/2 introduced by the confinement is 
exactly half of the number of colloids. The total interaction energy in the dimensionless 
form is: 

N N e -Kr tj N 

E = E E 77 - A E cos(27ra;;/0 + cos^vn/,//) , (1) 

i=l j^i r ijl l i=l 

where (xi,yi) are the Cartesian coordinates of colloid i in the plane, the distance 
between colloids i and j, \ j k is the Debye screening length. The unit of length I is the 
inter trap distance (see Fig. []~J, and the parameter A is a measure for the relative 
strength of the confinement over the electrostatic forces. The physical parameters 
governing the phase behaviour of the system are nl, A and the temperature T. At 
low temperatures and strong confinement each trap is occupied by exactly 2 colloids 
forming an elongated colloidal dimer. Due to the dimer-dimer electrostatic interactions, 
the dimers order into orient at ionally ordered colloidal molecular crystals. Upon raising 
the temperature first the orientational order of the dimers decays and eventually it 
becomes possible for single colloids to hop between the traps thus destroying the dimers 
and forming a modulated liquid state. 

We used parallel tempering Monte Carlo (PTMC) simulations to simulate M replicas 
of the system, each in the canonical ensemble, at temperatures Tj ranging from 
T\ = 10~ 4 to Tm = 1 allowing for configurational swaps between systems with adjacent 
temperatures [IH] (we chose a geometric temperature profile, thus Tj/T i+ i = const). 
Apart from single colloid moves and configurational swaps we also implemented dimer 
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rotations that become important at low temperatures. The colloids are treated as 
point-like charges and the periodic boundary conditions are applied to the simulation 
box. In Figj2] we show the energy histograms P(E) for our model system at lowest 
five temperatures and three system sizes. Overlap of P(E) between adjacent replicas 
at different temperatures allows for acceptance of the configuration swaps; in our 
simulations the target replica exchange probability is set to 0.2. At a given number of 
colloids N, the number of replicas is chosen accordingly. We found an empiric relation 
M ~ VoW and studied systems with (N, M) = (8, 20), (32, 40), (72, 60), (128, 80), 
and (200,100). From the upper inset in Figj2] we see that the width of P(E) (that 
is proportional to the standard deviation o E of the distribution) decreases as 1/y/N so 
that number of replicas M should indeed scale as \/N to maintain a constant overlap. 
We note that Tj (except 7\ and Tm) have different values for different system sizes N . 
We first examined the ground state phase diagram of dimers (n = 2) on a square lattice 
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Figure 2. a) Energy histograms (Gaussian distributions) P{E) calculated at five 
lowest temperatures T±...T^ in equilibrated PTMC simulations for three system sizes 
N = 32, 72, 1 28 (shifted for clarity). The values of the physical parameters are nl = 10 
and A — 0.6. Insets: Standard deviation erg of the distributions P{E) as a function of 
b) the number of colloids N at the lowest temperature T — T\ and c) the temperature 
T at the values of N displayed in the main plot. 

by PTMC simulations. The parameter space (nl, A) we have explored is depicted on 
Figj3] where several regions are highlighted. In most of the parameter space on the figure 
(right hand side) the ground state configuration of the colloids is antiferro-like (AF), 
meaning that the dimers directions in neighbouring traps are perpendicular. Structures 
different from AF are observed at low A and kI: typically, at vanishing confinement 
strength A, the colloids arrange into a hexagonal lattice and upon increasing A, the 
structure evolves via I-phases (checkerboard arrangement with neighbouring dimers' 
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relative angle less than tt/2, see [T4"IIT2] ) into AF arrangement (where the angle between 
neighbouring dimers is tt/2). At low values of the screening parameter kI, the electric 
double layers are of comparable size to the lattice spacing and many-body effects among 
colloids are likely to play an important role in orientational ordering. We have roughly 
estimated and highlighted the many-body dominated region with dark grey shading 
in FigJH Most of the non-AF phases we discussed above (including those reported in 
literature [10] where the parameters were kI = 2 and < A < 0.64) are in or close to 
this region, which means that their stability would need to be re-assessed by applying 
the many-body Poisson-Boltzmann theory. The behaviour in this regime is definitely 
interesting, however it is experimentally difficult to access and was not the focus of 
our present study. We would also like to stress that in a large part of the parameter 
space where the confinement is relatively weak the size of the "molecules" in the traps 
is significant with respect to the lattice spacing and the notion of molecule becomes 
obscure. We have provisionally marked this region with a light grey shading (please 
note that there is a certain degree of arbitrariness in defining the boundaries of such 
regions). To study the phase behaviour, we concentrated on the remaining part of the 
parameter space, where molecules are well defined and many-body effects are negligible. 
Such regime is also where most of the experimental studies [8, 9J have been performed. 
In this parameter regime the ground state at T = is AF phase. We have analysed the 
phase behaviour as a function of temperature at the points described by three sets of 
parameters (depicted by the diamond symbols in Fig^]): 0.1 < A < 10; nl = 10,6 and 
A = 1; 2 < kI < 10. 

3. Results 

In FigJH we show typical snapshot configurations at different temperatures across 
the melting transition. The right-most column shows the energy and moments of 
positional and angular distributions as a function of temperature. In the snapshots, 
one can observe AF configurations with long-range positional and orientational order, 
positionally ordered (exactly 2 colloids per trap) but orientationally disordered states 
and liquid-like disordered configurations. In order to quantify the positional order, 
we show normalized distributions of the trap occupancy P(N t ), which is narrow and 
peaked at N t = 2 at low enough temperatures, while it broadens around the average 
N t = 2 at higher temperatures. From these distributions it seems natural to extract 
a T-dependent order parameter as P{2), measuring the ratio of traps in the system 
with exactly 2 colloids. Any deviations from the saturated value P{2) = 1 signals 
the presence of positional melting. Smooth decrease of P{2) with increasing T is 
a sign of a crossover transition rather than a phase transition; this is additionally 
supported by the lack of symmetry breaking at the transition. In this respect, a 
crossover transition temperature T c \ cannot be rigorously defined, yet we decide to 
mark with T cl the beginning of crossover at the point where P(2) starts to deviate 
from 1. At that temperature on average the N tr — 2 traps still host exactly 2 colloids, 



Melting of colloidal molecular crystals 



6 



10 



A 



0.1 



0.01 



0.001 



1 

V 


1 i 1 


A 

V - 


A 

\ 




V - 


o 




O 


/\ /\ /\ 


<? 


<? - 


o 




z\ - 

<✓ : 


: \ ""■•••••■<> 




O 


r \ O 




On 




AF phase 
















: many-body\ ^. 


no molecules 




non\AF phases \ 

i i ,ii 


i 


i 



6 
Kl 



10 



Figure 3. T=0 phase diagram of dimers on a square lattice of traps. Line separating 
the AF region from the non-AF region is schematic and was calculated on a N = 32 
system using PTMC method. The diamond symbols are placed at the positions in 
the parameter space where the finite temperature behaviour has been explored. The 
regions of the parameter space where many-body interactions are dominant and where 
the molecules are not well defined are marked with the shades of grey. Note that phase 
diagram where cosine confinement is replaced by the parabolic trapping [14j is quite 
richer. 



while two of the traps exchange one colloid to host 1 and 3 colloids, respectively and 
the threshold value of the order parameter is P(2) = 1 — 2/N tr . From the results at 
nl = 10 and A = 0.6 we estimate T& ~ 2 x 10~ 2 . The so-defined quantity is showing 
negligible finite size effects. At low temperatures the colloidal dimers are orient ationally 
correlated. This is indirectly demonstrated in the bottom row of Fig. H] where the 
angle distributions P(<f>) of dimers at various temperatures [20] are depicted: at the 
lowest temperature two narrow peaks of equal height imply that half of the dimers lie 
horizontally (0 = 0) and the other half vertically (0 = tt/2), as expected for the AF 
phase. To characterize the above distributions with a meaningful scalar measure, we 
define the spread o v of the distribution P(<p): o~ 2 P = l/NunJ2i{Pi ~ 1) 2 ; where 
is the number of bins used to sample the variable (p. In order to have a quantity 
bounded by 1, we define the spreading parameter ap/o-p, where a P = ^/iV bin /2 — 1 is 
the zero temperature limit of ap in the AF phase. At large T the value of crp/crp is 
zero indicating orientationally disordered configuration and it starts increasing upon 
lowering the temperature below T c \. However, the spreading parameter is not a suitable 
order parameter for studying the details of the orientational melting transition, since it 
only deals with one-dimer statistics and therefore does not capture dimer-dimer angular 
correlations. To properly account for the angular correlations in dimeric coverings, 
we define a new order parameter analogue to the Fourier transform of the spin-spin 
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Figure 4. Results of the PTMC simulations of 72 colloids subject to periodic array 
of 36 traps at nl = 10 and A = 0.6. Top row: Typical colloidal configurations from 
the simulations at four selected temperatures. The right column shows the energy 
E(T) versus temperature T (the four vertical lines denote the temperatures of the 
configurations on left). Mid row: Distribution of the number of colloids per traps 
for the same parameters as in the top row. The right column shows the proportion 
of double-occupied traps P(2) as a function of temperature T. Bottom row: Angle 
distributions P(<fi) of dimers (or higher n-mers if there are more than two colloids 
per trap [20]) calculated at the same four temperatures as above. The right column 
shows the temperature dependence of the spreading parameter as a measure of angular 
anisotropy. 



correlator [21] . We start from the averaged static correlation function that measures 
average angle-angle correlation between two traps separated by 8 

1 iVtt 

o(«) = F E(l««(^-^> ( 2 ) 

iV *r i=l 

where N tT denotes the number of traps in a system and 0j measures orientation of the 
dimer in z-th trap and (...) denotes thermal averaging. Since dimers are not directed 
and therefore <j) and 0±7r indicate the same orientation, we use the absolute value in Eq. 
(J2]). The function 0(8) is bounded on [0, 1]; in order to compute its Fourier transform, 
we first expand its range to [—1, 1] by transforming it: 0(8) — > 20(8) — 1. The Fourier 
transform then reads: 

0(q) = ^£(20(<*)-l)e iq<5 (3) 

q being a 2D ordering wave vector. The temperature variation of O(q) across the 
melting transition from AF to modulated liquid is depicted in Fig|5l In an ideally 
ordered AF phase at T = 0, characterized by modulation q = qo, the order parameter 
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saturates, 0(q = (*7r,7r)) = 1. At high temperatures, on the contrary, any pair 
of traps is uncorrelated, which gives a vanishing O(q) for all non-zero wavevectors: 
0(q 7^ 0) oc 1/iVtr [22]- In Figj6] we compare the orientational and positional order- 




1.5 ^o.o 



Figure 5. The orientational order parameter O(q) as a function of the two-dimensional 
wave vector q at three temperatures across the melting transition (nl = 10, A = 1). At 
the lowest temperature (AF phase) the peak O(qo = (k,it)) reflects the modulation 
of the angular orientation of the dimers. Upon raising the temperature, this peak 
gradually decays and at high temperature (modulated liquid phase), where the dimers 
are uncorrelated, a peak at q = emerges. 



disorder transitions at kI = 10 for several confinement strengths A at three system 
sizes. In the upper row we plot the maximum value of the trap occupancy distribution 
P(2) as an order parameter for the positional order. We see that the melting transition is 
gradual, depends on the confinement strength and shows very little size dependence. In 
the lower row we show the temperature dependence of the orientational order parameter 
O(q ) (Eqj3]) using the antiferro ordering wave vector q = (tt, iy). It varies from 1 (AF 
phase) at low temperatures to (orientational disorder) at high T. In contrast to the 
positional melting, the transition is now sharp and it grows sharper with increasing 
system size. This behaviour is a signature of a second order phase transition with 
diverging correlation length. We also note that at the transition the system reduces 
translational symmetry with respect to the periodicity of the underlying trap potential. 
In FigfTl the same quantities are plotted for nl = 6. Very similar results have been 
obtained at all other parameter values kI and A marked on Figj3] (results not shown). 
We therefore conclude that there are -in agreement with previous results \1Q\- two 
transitions going on in the system upon cooling from the liquid phase: at T cl »s 10~ 2 a 
gradual crossover ordering transition to the positionally ordered structures sets in, and 
at T C 2 ~ 10~ 3 a second order phase transition from orientationally disordered to ordered 
state. In order to further strengthen the evidence for a second order phase transition, 
we have also calculated the specific heat as a function of temperature. The specific heat 
of a finite canonical system is given by the energy fluctuations: 

(E 2 ) - (£> 2 



C 



2^2 



(4) 



We have evaluated this expression and observed a narrow peak in C at T C 2 showing the 
expected system size dependence for a second order phase transition (in an infinite 
system the specific heat would diverge at T c2 ). In contrast, there is a broad size- 
independent peak positioned at T cl . In Fig. [SI we present the specific heat curves as a 
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Figure 6. Melting transitions of AF phase for kI — 10 and various trap strengths 
A = 0.2...1 and three system sizes. Upper row: the positional order parameter P{2) 
as a function of T. Lower row: the temperature dependence of the orientational 
order parameter O(qo = (tt, i")) defined in Eq. ([3]). 

N = 36 x 2 64x 2 100 x 2 




Figure 7. Same as Fig[6]for kI — 6 and trap strengths A from 0.1 to 5. 



function of temperature and system size at two extreme values of A and demonstrate 
that a maximum increase of both order parameters coincides with peak positions in 
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Figure 8. Upper row: Temperature dependence and finite size scaling of the specific 
heat per colloid C/N at kI = 10 and two values of A Lower row: The corresponding 
variation of the orientational O(qo) (dashed lines) and positional P(2) (dotted lines) 
order parameters. 



transition T cl shifts towards higher temperatures with increasing trap strength A, while 
the position of the phase transition T C 2 shows a weaker dependence on A. From 
analysing the finite size scaling (see the Supplementary material) of the two transition 
temperatures we conclude that both, T c \ and T c2 are finite in the thermodynamic limit 



4. Discussion 



To summarize, we present (Figf5J) phase diagrams showing both transition temperatures 
versus kI at fixed confinement strength and versus A at fixed kI. The crossover transition 
T c i has been somewhat arbitrary defined as the temperature at which the positional 
melting starts and the value of the order parameter P(2) drops from 1 to 1 — 2/N tz . 
On the other hand, T c2 was extracted as the point of the steepest decrease of the 
orientational order parameter O(qo) exactly matching the position of the diverging 
peak in the specific heat C (FigJBJ). With increasing screening strength kI, both critical 
temperatures decrease (see left panel in Fig. [9]). This is due to weaker inter-trap 
interactions that directly control orientational ordering (T c2 ) and also, in addition to 
trapping forces, average distances between the dimers, which affects positional ordering 
(T c i). Increasing A has a similar effect on T c2 (see middle and right panels in Fig. 
[9} since larger A means smaller dimers and hence weaker inter-trap forces. On the 
contrary, positional order is of course stronger at larger A, therefore higher temperatures 
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are needed to destroy it. At large enough values of kI and A we observe the two-stage 
melting with the phase transition from ordered solid to partially ordered solid (without 
orientational order) and the crossover melting of the positional order from partially 
ordered solid to modulated liquid. Interestingly, at longer screening lengths (smaller 
kI) and weaker confinement the crossover transition disappears and the ordered solid 
melts directly into a modulated liquid via a second order phase transition. This can 
also be observed in Fig. [7] (first row, black line) where the P(2) curve exhibits a sudden 
change in the slope, precisely at T c2 . The phase diagram reported here is qualitatively 

A=l „ Kl=6 Kl=10 
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io- 2 

T 

IO" 3 



2 4 6 8 10 1U 0.1 1 10 0.1 1 10 1U 
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Figure 9. Phase behaviour of dimers on the square lattice as calculated in the 
largest system with A?~ = 200 colloids: a) Transition temperatures versus kI at fixed 
confinement strength A = 1; b) Transition temperatures versus A at fixed screening 
kI = 6 and c) kI = 10. 




consistent with the result of Langevin simulations reported in pU]. The screening 
lengths considered in our work (nl = 6,10) are much shorter than in [10] (kI = 2), 
therefore the phase diagrams in Figj9] are shifted to lower values of the confinement 
A as compared to Fig.4 in |1Q| . In [13] spin models were constructed to analytically 
describe orientational ordering of colloidal molecules on periodic lattices. Remarkably, 
their analytically obtained orientational melting transition for trimers on triangular 
lattice coincides with the corresponding transition for dimers on square lattices in 
[10J. Crystaline structures of dimers on triangular lattice were also examined in [T3] 
showing interesting transitions between herringbone, ferromagnetic, "Japanese 6 in 1", 
paramagnetic and antiferromagnetic ordering. However, several approximations had to 
be made in order to reduce the colloidal system to a Potts-like model: they assumed 
rigid dimers placed on the lattice points interacting with nearest neighbours only and 
restricted to discrete orientations compatible with the lattice symmetry. In |14| we 
have already discussed that the discrete angle approximation does not apply to the 
experimental system of colloids trapped by laser tweezers. Further to that, while 
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the nearest neighbour approximation seems well justified for triangular lattices, it is 
qualitatively wrong on square substrates [H], therefore it seems difficult to extend the 
conclusions of [13] to colloidal systems and to square-like patterned substrates. 
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exponent more accurately. The obtained value lies in between the spin-model values for Ising 
model (a = 0) and 3-state Potts model (a = 1/3) [24]. The corresponding exponent for the 
order parameter O(qo) proved to be more difficult to asses by finite size scaling. Roughly, the 
value seems to be (3 « 0.2 in contrast to the spin model values (Ising: 1/8 and 3-state Potts: 
1/9). The more elaborate evaluation of the critical exponents is out of the scope of the present 
manuscript. 

[24] F. Y. Wu, Reviews of Modern Physics, 54, 235 (1982). 



